## Table and figures script file
## Time-stamp: <2016-9-13 TH>
## Contact: hasegaw.tomoko@nies.go.jp

getwd()
setwd(getwd())
#setwd("D:\thasegawa\ExtremeEvent\R")
library(reshape2)
library(ggplot2)
library(plyr)
library(grid)
library(scales)

colnames <- c("VAR","Year","Region","SSP","GCM","RCP","CO2fert","SYR","N","Crop","Value")
yearname <- c("2050")

fsize <- 15
MyThemeLine <- theme_grey() +
  theme(
    panel.border=element_rect(fill=NA),
    title=element_text(size=15,face="bold"),
    axis.title=element_text(size=15,face="bold"),
    axis.text = element_text(hjust=1,size = 15, angle = 0),
    axis.line=element_line(colour="black"),
    panel.background=element_rect(fill = "white"),
    #    panel.background=element_blank(),
    #    panel.grid.major=element_line(linetype="dashed",colour="grey",size=0.5),
    panel.grid.major=element_blank(),
    strip.background=element_rect(fill="white", colour="white"),
    strip.text.x = element_text(size=15, colour = "black", angle = 0,face="bold"),
    legend.title = element_text(size=15),
    legend.text = element_text(size=13)
  )

factorlabel <- function(y){
  levels(y$VAR)[levels(y$VAR)=="FAREA"] <- "Area"
  levels(y$VAR)[levels(y$VAR)=="FFEED"] <- "Feed"
  levels(y$VAR)[levels(y$VAR)=="FNETT"] <- "Trade"
  levels(y$VAR)[levels(y$VAR)=="FOTH"] <- "OthCons"
  levels(y$VAR)[levels(y$VAR)=="FYILD"] <- "Yield"
  levels(y$VAR)[levels(y$VAR)=="error"] <- "error"
  z <- y
  return(z)
}


pdata <- read.table("../data/Data_decomp.txt")
names(pdata) <- colnames

pdata$RCP <- factor(pdata$RCP, levels=c("RCP85","RCP26"))
pdata$VAR <- factor(pdata$VAR, levels=c("dFOOD","FAREA","FYILD","FNETT","FFEED","FOTH","error"))

dirname <- paste('../output/decomp', sep = '')
dir.create(dirname)

labelname <- c("Change in calorie intake [kcal/person/day]")

y <- pdata[pdata$Year==yearname & pdata$Region=="IND"  & pdata$Crop=="WHT" & pdata$SSP=="SSP2", c(-2,-4)]

names(y) <- c("VAR","Region","GCM","RCP","CO2fert","SYR","N","Crop","Value")  #Rename the column
interact<-interaction(y$RCP, y$CO2fert, sep=" : ")

minv <- min(y[9])
maxv <- max(y[9])

y1 <- ddply(y,c("VAR","interact"),transform,v0=min(Value))
y2 <- ddply(y1,c("VAR","interact"),transform,v25=quantile(Value,0.25))
y3 <- ddply(y2,c("VAR","interact"),transform,v50=median(Value))
y4 <- ddply(y3,c("VAR","interact"),transform,v75=quantile(Value,0.75))
y5 <- ddply(y4,c("VAR","interact"),transform,v100=max(Value))

z <- factorlabel(y5)

g <- ggplot() + 
  geom_boxplot(data=z, aes(x=interact,y=Value, color=VAR, ymin = v0, lower = v25, middle = v50, upper = v75, ymax = v100), fill="grey80", stat="identity") +
  ylim(minv-abs(minv)*0.01, maxv+abs(maxv)*0.01) +
  #	ggtitle(mainlabel)
  ylab(labelname) + xlab("") +	MyThemeLine

plot(g)

name <- paste('../output/decomp/',yearname,'_cp4.png', sep = '')

ggsave(g,filename=name, dpi=250,width=10, height=6)
